读取ROOT的tree数据，进行逐事件分析
一、学习目的¶
1、学习进行逐事件分析的方法 2、掌握中子探测器位置刻度、飞行时间刻度的方法 3、从实验测量值 tu,td,Qu,Qd 推出 学习数据分析的思路
二、具体的代码实施过程-李老师¶
In [1]:
TCanvas *c1=new TCanvas();//* 在ROOT环境下可省略
TH1D *hTOF= new TH1D("hTOF","Time of flight", 1000,0,100);
打开文件，得到TTree的指针¶
In [2]:
TFile *Input_File=new TFile("tree.root");//打开ROOT文件
if (Input_File->IsZombie()) {
cout << "Error Open File" << endl;
exit(-1);
}
Input_File->cd();
TTree *tree=(TTree*)Input_File->Get("tree");//得到名字为“tree”的TTree指针
In [3]:
//声明所需要的参数，即需要输出的参数
//* 入射粒子原始信息：粒子种类、能量，粒子在探测器上的入射位置、入射深度等
//* 探测器信息: 两端时间信息 tu 、 td ，能量沉积信息 Qu , Qd
Double_t X,Energy; //应该是均匀分布
Int_t PID; //pid=1 is neutron、pid=0 is gamma
Double_t T_u,T_d,Q_u,Q_d;
Double_t TOF,cTOF; //TOF:粒子实际飞行时间，cTOF：计算得到的TOF
声明tree的Branch变量¶
In [4]:
//3. 将变量指向对应Branch的地址
tree->SetBranchAddress("cTOF",&cTOF);//将ROOT文件内tree内名为"ctof"的branch的数据的指针指向ctof的变量。
tree->SetBranchAddress("TOF",&TOF);
tree->SetBranchAddress("PID",&PID);
tree->SetBranchAddress("T_u",&T_u);
tree->SetBranchAddress("T_d",&T_d);
tree->SetBranchAddress("Q_u",&Q_u);
tree->SetBranchAddress("Q_d",&Q_d);
为将新数据写入新的ROOT文件做准备¶
In [5]:
Double_t T2_u,T2_d;
Double_t cTOF2,cX;
TFile *Output_File2=new TFile("tree2.root","recreate");
///Output_Tree->cd();
TTree *Output_Tree2=new TTree("tree2","tree2 structure");
Output_Tree2->Branch("T2_u",&T2_u,"T2_u/D");
Output_Tree2->Branch("T2_d",&T2_d,"T2_d/D");
Output_Tree2->Branch("cX",&cX,"cX/D");
Output_Tree2->Branch("cTOF2",&cTOF2,"cTOF2/D");
逐事件读取tree的branch数据¶
In [6]:
Long64_t nentries=tree->GetEntries();//得到tree的事件总数
for(Long64_t jentry=0; jentry<nentries; jentry++) {//对事件进行遍历
tree->GetEntry(jentry);//将第jentry个事件数据填入对应变量(步骤3.中指向的变量)，每次变量值会变成当前事件对应的数据。
hTOF->Fill(cTOF);
// calculate new parameters
T2_u=T_u;
T2_d=T_d;
cTOF2=(T2_u+T2_d)/2;
cX=(T2_u-T2_d)/2;
Output_Tree2->Fill();//fill new parameter to TTree* opt
}
hTOF->Draw();
Input_File->Close();
Output_Tree2->Write();
Output_File2->Close();
c1->Draw();
三、验证是否存储正确¶
In [7]:
%jsroot on
TFile *Input2_File=new TFile("tree2.root");
Input2_File->ls()
In [8]:
tree2->Print()
In [9]:
c1->Clear();//* 在ROOT环境下可省略
tree2->Draw("T2_d-T2_u>>htx(500,-20,50)");//位置一维分布
c1->Draw();//* 在ROOT环境下可省略
In [10]:
c1->Clear();
tree2->Draw("cTOF2>>hh(1000,0,100)");//实际飞行时间
c1->Draw();
In [11]:
c1->Clear();
tree2->Draw("cX>>hh(1000,-50,50)");//实际飞行时间
c1->Draw();
In [12]:
c1->Clear();
tree2->Draw("T2_u-T2_d:cX");//实际飞行时间
c1->Draw();
四、具体的代码实施过程-吴师兄¶
In [13]:
TFile *Input_File3=new TFile("tree.root");//打开ROOT文件
if (Input_File3->IsZombie()) {
cout << "Error Open File" << endl;
exit(-1);
}
Input_File->cd();
TTree *tree3=(TTree*)Input_File3->Get("tree");//得到名字为“tree”的TTree指针
tree3->SetBranchAddress("cTOF",&cTOF);//将ROOT文件内tree内名为"ctof"的branch的数据的指针指向ctof的变量。
tree3->SetBranchAddress("TOF",&TOF);
tree3->SetBranchAddress("PID",&PID);
tree3->SetBranchAddress("T_u",&T_u);
tree3->SetBranchAddress("T_d",&T_d);
tree3->SetBranchAddress("Q_u",&Q_u);
tree3->SetBranchAddress("Q_d",&Q_d);
tree3->SetBranchAddress("X",&X);
In [14]:
Double_t T3_u,T3_d,Q3_u,Q3_d;
TFile *Output_File3=new TFile("tree3.root","recreate");
///Output_Tree->cd();
TTree *Output_Tree3=new TTree("tree3","tree3 structure");
Output_Tree3->Branch("T3_u",&T3_u,"T3_u/D");
Output_Tree3->Branch("T3_d",&T3_d,"T3_d/D");
Output_Tree3->Branch("Q3_u",&Q3_u,"Q3_u/D");
Output_Tree3->Branch("Q3_d",&Q3_d,"Q3_d/D");
TH1D *tdiff=new TH1D("tdiff","T3_d-T3_u",350,-20,50);
In [15]:
c1->Clear();
for(Long64_t jentry=0; jentry<nentries; jentry++) {//对事件进行遍历
tree3->GetEntry(jentry);//将第jentry个事件数据填入对应变量(步骤3.中指向的变量)，每次变量值会变成当前事件对应的数据。
// calculate new parameters
T3_u=T_u;
T3_d=T_d;
Q3_u=Q_u;
Q3_d=Q_d;
cX=(T3_d-T3_u)/2;
tdiff->Fill(cX*2);
}
tdiff->Draw();
c1->Draw();
In [16]:
TH1D *dtd=new TH1D("dt_d","dt/dx",141,-20.25,50.25);
In [17]:
for(int i=1;i<=tdiff->GetNbinsX();i++) {
Double_t df=tdiff->GetBinContent(i+1)-tdiff->GetBinContent(i); //在bin中循环，在bin中得到内容
dt_d->Fill(tdiff->GetBinLowEdge(i+1),df); //在root里面自带误差传递
}
//dt_d->Sumw2(0);//不显示传递误差 =等价于在Draw中加“hist”
dt_d->Draw("hist");
dt_d->Fit("gaus","","",-14,-10.5);//txl，需要明白，边界
c1->Draw();
拟合存在问题的原因是变化趋势的缘故,gaus函数不能拟合负数？
In [18]:
dt_d->Fit("gaus","","",39,45);//为什么拟合有问题？
c1->Draw();
In [19]:
TF1 *f1 = new TF1("f1","[0]*TMath::Exp(-0.5*((x-[1])/[2])^2)",40,43);//定义函数的方法 TF1，定义函数的方法
//需要按照顺序来写
//gaus:f(x) = p0*exp(-0.5*((x-p1)/p2)^2)
In [20]:
//进行参数拟合时，设置合理的初始参数至关重要！
c1->SetLogy(0);
f1->SetParameter(0,-350); //大致估计范围
f1->SetParameter(1,41.5); //参数的调节不太懂
f1->SetParameter(2,0.5);
dt_d->Fit("f1","R");
dt_d->Draw();
c1->Draw();
In [21]:
const Double_t T3R=41.5834,T3L= -11.7698;
In [22]:
const Double_t Dis=500.,Len=100.,del_Dis=5.; //单位为cm，
Double_t t_xoff=(T3R+T3L)/2; //计算度越时间，即传输时间
Double_t t_xon=(T3R-T3L)/2;
In [23]:
//位置刻度
In [24]:
TH1D *htx=new TH1D("htx","htx",500,-120,120);
TH2D *hdx=new TH2D("hdx","htx-hx:hx",100,-20,20,500,-120,120);
刻度后的结果¶
In [25]:
for(Long64_t jentry=0; jentry<nentries; jentry++) {//对事件进行遍历
tree3->GetEntry(jentry);
T3_u=T_u;
T3_d=T_d;
Q3_u=Q_u;
Q3_d=Q_d;
Double_t tx=Len/t_xon*(T3_d-T3_u-t_xoff);
htx->Fill(tx);
hdx->Fill(tx-X,X);//difference 检查依赖关系，要求没有关系
}
htx->Draw();
c1->Draw();
In [26]:
hdx->Draw("colz");//为一条竖线，没有关联
c1->Draw();
In [27]:
TH1D *hdx1=hdx->ProjectionX("projx of hdx");//投影
hdx1->Draw();
hdx1->Fit("gaus");//高斯拟合
c1->Draw();
In [28]:
c1->Clear();
tree3->Draw("cTOF>>(1000,20,120)");
c1->Draw();
对TOF进行飞行距离修正¶
In [29]:
//以上结果没有问题
TH2D *hgtofx=new TH2D("hgtofx","hgtofx",100,-120,120,100,7,9);
TH1D *hgctof=new TH1D("hgctof","hgctof",100,39,48);
In [30]:
for(Long64_t jentry=0; jentry<nentries; jentry++) {//对事件进行遍历
tree3->GetEntry(jentry);
T3_u=T_u;
T3_d=T_d;
Q3_u=Q_u;
Q3_d=Q_d;
Double_t tx=Len/t_xon*(T3_d-T3_u-t_xoff);
if(cTOF>40&& cTOF<44) { //选择gamma
Double_t d=TMath::Sqrt(502.5*502.5+tx*tx);
Double_t cTOFa=(cTOF)/d*100.;//normalized to 500cm
// cout<<cTOFa<<endl;
hgtofx->Fill(tx,cTOFa);
if(abs(tx)<5) hgctof->Fill(cTOF);//gamma hits the center of the det.
}
}
hgtofx->Draw("colz");
c1->Draw();
In [31]:
//击中中心
c1->Clear();
hgctof->Draw();
c1->Draw();
In [32]:
TH1D *hC=new TH1D("hC","",300,-40,-10);
TH2D *hCx=new TH2D("hCx","",100,-105,105,300,-26,-23);
In [33]:
for(Long64_t jentry=0; jentry<nentries; jentry++) {//对事件进行遍历
tree3->GetEntry(jentry);
T3_u=T_u;
T3_d=T_d;
Q3_u=Q_u;
Q3_d=Q_d;
Double_t tx=Len/t_xon*(T3_d-T3_u-t_xoff);
if(tx>-100 && tx<100 && cTOF>40&& cTOF<44) {//gamma
Double_t d=TMath::Sqrt(502.5*502.5+tx*tx);
hC->Fill(3.333*d*0.01-cTOF);
hCx->Fill(tx,3.333*d*0.01-cTOF);
}
}
hC->Draw("");
hC->Fit("gaus","","",-25,-24);
c1->Draw();
In [34]:
Double_t C=-24.4947;
In [35]:
hCx->Draw("colz");//检验C与x是否有关联
c1->Draw();
In [36]:
TH1D *htof=new TH1D("tof","",500,0,30);//real TOF, normalized to 500cm
TH1D *htofc=new TH1D("tofc","",500,0,30);//calculated TOF, normalized to 500cm
TH2D *htofcx=new TH2D("htofcx","htofcx",100,-120,120,100,0,5);
In [37]:
for(Long64_t jentry=0; jentry<nentries; jentry++) {//对事件进行遍历
tree3->GetEntry(jentry);
T3_u=T_u;
T3_d=T_d;
Q3_u=Q_u;
Q3_d=Q_d;
Double_t tx=Len/t_xon*(T3_d-T3_u-t_xoff);
Double_t d=TMath::Sqrt(502.5*502.5+tx*tx);
Double_t TOF_d100=TOF/d*100;
Double_t cTOF_d100=(cTOF+C)/d*100;
// cout<<TOF_d100<<endl;
htof->Fill(TOF_d100);
htofc->Fill(cTOF_d100);
if(cTOF>40&& cTOF<42)
htofcx->Fill(tx,(cTOF+C)/d*100);
}
In [38]:
c1->SetLogy();
htofc->Draw();//修正后的飞行时间谱
htof->SetLineColor(kRed);
htof->Draw("same");
c1->Draw();
In [39]:
c1->SetLogy(0);
htofcx->Draw("colz");
c1->Draw();
In [40]:
TH1F *hce=new TH1F("hce","",500,0,200);//real TOF, normalized to 500cm
//TH1D *htofc=new TH1D("tofc","",500,0,30);//calculated TOF, normalized to 500cm
//TH2D *htofcx=new TH2D("htofcx","htofcx",100,-120,120,100,0,5);
In [41]:
for(Long64_t jentry=0; jentry<nentries; jentry++) {//对事件进行遍历
tree3->GetEntry(jentry);
T3_u=T_u;
T3_d=T_d;
Q3_u=Q_u;
Q3_d=Q_d;
Double_t tx=Len/t_xon*(T3_d-T3_u-t_xoff);
Double_t d=TMath::Sqrt(502.5*502.5+tx*tx);
Double_t TOF_d100=TOF/d*100;
Double_t cTOF_d100=(cTOF+C)/d*500;
Energy=72.29824/cTOF_d100*(d*0.01)*72.29824/cTOF_d100*(d*0.01);
hce->Fill(Energy);
}
In [42]:
hce->Draw();
hce->Fit("gaus","","",40,160);
c1->Draw();
In [43]:
!jupyter nbconvert 1.2_Read_tree.ipynb --to html
In [ ]: